1 Présentation du projet

Ce projet s’inscrit dans le cadre de la formation DSA, et évalue l’apprentissage de R réalisé auprès de Robin RYDER.

Le thème de ce projet est libre.

> J’ai choisi d’étudier une série d’articles qui m’ont frappé sur la comparaison des approches classiques par modèles linéaires généralisés GLM) et par machine learning (ci après ML) dans le cadre de la tarification en assurance non-vie et la manière dont elles pouvaient se compléter.

Dans le suite de ce rapport, on ne s’intéressera qu’à la modélisation de la fréquence.

Soit $N_{i}$ le nombre de sinistre de la police ${i}$, alors dans le modèle Poisson homogène on suppose $N_{i} \sim \mathscr{P}(\lambda \nu_{i})$ où $\nu_{i}$ est l’exposition au risque de la police $i$, $\lambda$ est la fréquence de sinistre annuelle (homogène sur le portefeuille)et $\mathscr{P}$ la loi de Poisson. On pose $Y_{i} = \frac{N_{i}}{\nu_{i}}$ la fréquence annuelle des sinistres pour la police $i$ et $ = (x_{1},x_{2}, , x_{q}) $ les $q$ caractéristiques de la police $i$.

En tarification, on suppose dans les faits que $\lambda$ dépende de $\mathbf{x}$, de sorte que $\mathbb{P}[Y=k/\nu] = \mathbb{P}[N=k] = e^{-\lambda(\mathbf{x})\nu} \frac{(\lambda(\mathbf{x})\nu)^{k}}{k!}$ avec $ln(\lambda(\mathbf{x})) = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + \cdots + \beta_{q}x_{q} \overset{def}{=} \langle \mathbf{\beta}, \mathbf{x} \rangle$ en complétant $\mathbf{x}$ par $x_{0} = 1$.

Lorsqu’on utilise les approches classiques de GLM, on estime $\hat{\lambda}$ en utilisant l’estimateur du maximum de vraisemblance $\hat{\beta}$ de $\beta$.

Cet estimateur peut être obtenu soit en maximisant la log vraisemblance, soit de manière équivalente en minimisant la déviance de Poisson.

Une première adaptation des modèles classiques consiste donc à définir la déviance de poisson comme fonction de perte des algorithmes de ML.

C’est notamment l’approche préconisée par dans Data Analytics for Non-Life Insurance Pricing de Mario V. Wuthrich et Christoph Buser

En complément, lorsqu’on utilise un GLM incluant un intercept et sa fonction de lien canonique, l’utilisatiuon d’un estimateur du maximum de vraisemblance apporte de bonnes propriétés pour son usage en tarification dont celle dite de propriété d’équilibre (Property 2.4 p33 de Data Analytics for Non-Life Insurance Pricing), qui précise que :

$\displaystyle\sum_{i=1}^{n} \nu_i\hat{\lambda}(\mathbf{x_{i}}) = \displaystyle\sum_{i=1}^{n} \nu_i exp\langle \mathbf{\hat{\beta}}, \mathbf{x_{i}} \rangle = \displaystyle\sum_{i=1}^{n} N_{i}$

De sorte que le nombre total de sinistres estimés est égale au nombre total de sinistres observés.

Cette propriété est même plus étendue dans certains modèles GLM, notamment lorsque le modèle poissonien avec lien log (exemple (b) p423 de A systematic relationship between minimum bias and gener- alized linear models) : dans ce cas l’équilibre existe au globaldu portefeuille, mais aussi pour chaque classe des variables catégorielles incluses dans le modèle.

Ces propriétés d’équilibres sont critiques pour un assureur dans un exercice de tarification, puisqu’elles assurent qu’une structure tarifaire basée sur cette approche lui permettra d’équilibrer (au moins théoriquement) ses encaissement avec les sinistres qu’il s’est engagé à indemniser. En outre, ces approches sont très intuitives et ont en partie contribué au succès des GLM en assurance.

Or il semble que l'utilisation directe de modèles de ML, mêmes adaptés pour minimiser la déviance, ne présentent pas ces garanties.

En effet, si les modèles de ML permettent une meilleurs corrélation avec la réponse individuelle de chaque police, aucune garantie n’est expressément introduite au global du portefeuille ce qui peut produire des divergences, potentiellement importantes, entre la somme des primes demandées et les sinistres que l’assureur s’est engagé à régler.

Cet effet est attribué par M Wütrich aux conditions d’arrêt précoce des algorithmes d’optimisation utilisés par les modèles ML

C’est l’objet de l’article Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning de Michel Denuit, Arthur Charpentier et Julien Trufin.

Cet article repart des écarts empiriquement constatés par M Wütrich et leur donne une

Data Analytics for Non-Life Insurance Pricing de Mario V. Wuthrich et Christoph Buser Cet article détaille les mathématiques derrières l’utilisation des modèles GLM et des principaux modèles ML, y compris les réseaux de neurones de type feed forward
Peeking into the black box de Christian Lorentzen et Michael Mayer Cet article met en avant une utilisation des techniques d’explicabilité des modèles ML pour comprendre la manière dont ces modèles combinent les informations en entrée pour arriver à leurs prédictions. Ce faisant le praticien est invité à :
- identifier le gain marginal du modèle ML par rapport au modèle GLM généralement en place et en cas d’écart significatif à
- identifier les éléments appris par le modèles ML (effets non linéaires adns les variables simples, principales intéractions…) qui peuvent être réintégrés dans le modè-le GLM afin d’en réduire l’écart en termes de performance
Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning de Michel Denuit, Arthur Charpentier et Julien Trufin Cet article est le plus intriguant des 3

nous avons eu l’occasion d’analyser un jeu de données de tarification MRH à l’occasion d’un Hackathon.

Du fait des contraintes propres au Hackathon je n’avais pas eu l’occasion de finaliser l’analyse de ce jeu de données. Par ailleurs, pour le hackathon, le choix s’était porté sur le logiciel Python.

Ce TD est donc l’occasion de poursuivre l’analyse de ce jeu de données. Toutefois, cet exercice est aussi l’occasion d’explorer les apports de deux articles :

Peeking into the black box de Christian Lorentzen et Michael Mayer
Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning de Michel Denuit, Arthur Charpentier et Julien Trufin

Ces articles s’appuient sur les travaux de M Wütrich qui proposaient d’adapter la fonction de perte des modèles de machine learning au cas poissonien (déviance de poisson) pour pouvoir appliquer les techniques de ML à des problèmes d’estimation de la fréquence de sinistres.

Cette approche diffère des consignes du Hackathon qui préconisaient de minimiser le RMSE. En effet, cette fonction d’erreur suppose une distribution normale, ce qui dans les cas de tarfication où la fréquence de sinistre est faible n’est pas réaliste.

Dans le suite de ce rapport, on ne s’intéressera qu’à la modélisation de la fréquence.

Lorsqu’on utilise les approches classiques de GLM, on estime $\hat{\lambda}$ en utilisant l’estimateur du maximum de vraisemblance $\hat{\beta}$ de $\beta$.

Cet estimateur peut être obtenu soit en maximisant la log vraisemblance, soit en minimisant la déviance de Poisson.

Une première adaptation des modèles classiques consiste donc à définir la déviance de poisson comme fonction de perte des algorithmes de ML.

En complément, lorsqu’on utilise un GLM incluant un intercept et sa fonction de lien canonique, l’utilisatiuon d’un estimateur du maximum de vraisemblance apporte de bonnes propriétés pour son usage en tarification dont celle dite d’équilibre (Property 2.4 p33 de Data Analytics for Non-Life Insurance Pricing), qui précise que :

$\displaystyle\sum_{i=1}^{n} \nu_i\hat{\lambda}(\mathbf{x_{i}}) = \displaystyle\sum_{i=1}^{n} \nu_i exp\langle \mathbf{\hat{\beta}}, \mathbf{x_{i}} \rangle = \displaystyle\sum_{i=1}^{n} N_{i}$

De sorte que le nombre total de sinistres estimés est égale au nombre total de sinistres observés.

Cette propriété semble même plus étendue dans les modèles GLM puisque cet équilibre se retrouve au global du portefeuille et pour toute modalité des variables catégorielles (dixit l’introduction de Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning)

Il semble toutefois que ce ne soit pas le cas pour les modèles de ML, même lorque ceux-ci sont paramétrer pour minimiser la déviance de Poisson. Wütrich attribue cet effet aux méthodes d’arrêt précoce des algorithme d’optimisation utilisés. L’explication qualitative produite est que les modèles de ML, par leur plus grande souplesse, permettent un meilleurajustement aux réponses individuelles. Toutefois, ceux-ci n’incluent pas de contraintes d’équilibre au niveau du portefeuille, de sorte que ce meilleur ajustement local peut créer un biais (potentiellement important) au niveau global.

Si cette propriété n’est pas respectée, il s’agit là d’un frein important à l’utilisation des modèles ML (non modifiés) en tarification, puisque l’assureur ne disposerait plus de garanties d’équilibre financier global entre son encaissement et les sinsitres qu’il s’est engagé à régler.

C’est le sens de l’article Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning) qui propose de se servir des primes estimées lors d’ue première étape par les modèles ML comme d’entrants pour une deuxième étape d’autocalibration locale. Cette méthode consiste à ajuster des modèles GLM locaux réduits à un intercept sur des ensembles de points qui sont définis à partir de la proximité des primes estimées par ML (l’étape de ML rassemble les observations qui se ressemblent, l’étape locale GLM lisse les écarts entre observations voisines pour retrouver les propriétés d’équilibre).

Soit $\mu_{i}(x_{i}) = \mathbb{E}[Y_{i}|\mathbf{X} = \mathbf{x_{i}}]$ Dans la

Ce TD est donc l’occasion de poursuivre l’analyse de ce jeu de données et d’explorer l’usage de packages et d’approches développés par MV Wütrich et disponible sur le site de l’association suisse des actuaires.

2 Analyse

2.1 Mise en place de l’environnement

Le code suivant est un ensemble d’utilitaire pour naviguer dans les répertoires relatifs du projet et installer ou charger les packages nécessaires à l’analyse :

if (!require("here")){ 
  install.packages("here") 
  library("here")
}

set_here

LoadPackage <- function (load.lib=c("")) {

  install.lib<-load.lib[!load.lib %in% installed.packages()]
  for(lib in install.lib) install.packages(lib,dependencies=TRUE)
  sapply(load.lib,require,character=TRUE)
  
}


LoadPackage(c('ggplot2', 'dplyr', 'formattable', 'DT', 'tidyr', 'caret', 'plotly', 'xgboost', 'flashlight', 'MetricsWeighted', 'corrplot', 'lsr', 'h2o'))

2.2 code

3 Présentation du jeu de données

Le jeu de données est constitué de 4 fichiers :

list.files(path=here('data', 'raw'))

## [1] "expo_test.csv"  "expo_train.csv" "primes2019.csv" "sin_train.csv"

Le projet suit un portefeuille d’assurance MRH suivi sur plusieurs années de 2016 à 2018 inclus. L’objectif initial du Hackathon est d’estimer les primes 2019 sur un échantillon de contrats.

3.1 expo_train

Le fichier expo_train.csv contient la liste des contrats en risques historiquement suivis :

expo_train = read.table(file = here('data', 'raw', 'expo_train.csv'), header=T, sep=',', dec='.', encoding = 'UTF-8', stringsAsFactors = F)

datatable(head(expo_train))

str(expo_train)

## 'data.frame':    155651 obs. of  19 variables:
##  $ X                  : int  384538 441079 119668 5124 130065 450830 151401 296526 71801 360445 ...
##  $ EXPO               : num  1 0.825 1 1 1 ...
##  $ FORMULE            : chr  "MEDIUM" "CONFORT" "ESSENTIEL" "ESSENTIEL" ...
##  $ TYPE_RESIDENCE     : chr  "PRINCIPALE" "PRINCIPALE" "PRINCIPALE" "SECONDAIRE" ...
##  $ TYPE_HABITATION    : chr  "APPARTEMENT" "MAISON" "APPARTEMENT" "MAISON" ...
##  $ NB_PIECES          : int  1 NA 3 2 1 2 2 2 1 3 ...
##  $ SITUATION_JURIDIQUE: chr  "PROPRIO" "PROPRIO" "LOCATAIRE" "LOCATAIRE" ...
##  $ NIVEAU_JURIDIQUE   : chr  "JUR1" "JUR1" "JUR1" "JUR1" ...
##  $ VALEUR_DES_BIENS   : num  3500 0 35000 9000 20000 9000 20000 9000 3500 0 ...
##  $ OBJETS_DE_VALEUR   : chr  "NIVEAU_1" "NIVEAU_1" "NIVEAU_1" "NIVEAU_1" ...
##  $ ZONIER             : chr  "B40" "A11" "B32" "C24" ...
##  $ NBSIN_TYPE1_AN1    : int  0 0 0 0 0 0 1 0 0 0 ...
##  $ NBSIN_TYPE1_AN2    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ NBSIN_TYPE1_AN3    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ NBSIN_TYPE2_AN1    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ NBSIN_TYPE2_AN2    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ NBSIN_TYPE2_AN3    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ id                 : int  5 9 11 13 14 15 20 21 28 30 ...
##  $ ANNEE              : int  2017 2018 2017 2018 2017 2017 2018 2016 2018 2016 ...

Les données ont été prétraitées pour disposer d’une ligne par période de risque annuelle homogène.

Pour des raisons d’anonymisation, les identifiants contrats ont été supprimés de sorte qu’il n’est pas possible de suivre l’évolution du risque d’une période sur l’autre. Pour la même raison, il ne sera pas possible dans la construction des jeux de validation et de test de s’assurer que toute l’expérience d’un même assuré est bien intégralement dans un jeu d’entraînement, de validation ou de test.

La signification des colonnes présentes est la suivante :

3.1.1 `X`:

variable non nommée. Il s’agit d’un artefact du processus de création du fichier initial. On ne la prend pas en compte

3.1.2 `EXPO` :

exposition en année risque du contrat. Sa valeur précalculée pour chaque ligne est comprise entre 0 et 1

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.002732 0.811475 1.000000 0.837278 1.000000 1.000000

3.1.3 `FORMULE` :

variable catégorielle codant la formule du contrat comprenant 3 niveaux MEDIUM, ESSENTIEL et CONFORT

FORMULE	nb_lignes	EXPO	pct_EXPO
ALL_INCLUDE	6,671.00	4,759.51	3.65%
CONFORT	87,603.00	73,930.57	56.73%
ESSENTIEL	22,364.00	19,546.44	15.00%
MEDIUM	39,013.00	32,086.59	24.62%

3.1.4 `TYPE_RESIDENCE` :

variable catégorielle codant le fait que le bien est une résidence PRINCIPALE, ou SECONDAIRE

TYPE_RESIDENCE	nb_lignes	EXPO	pct_EXPO
PRINCIPALE	130,727.00	112,793.55	86.55%
SECONDAIRE	24,924.00	17,529.57	13.45%

3.1.5 `TYPE_HABITATION` :

variable catégorielle codant le fait que le bien est un APPARTEMENT, ou une MAISON

TYPE_HABITATION	nb_lignes	EXPO	pct_EXPO
APPARTEMENT	71,123.00	61,562.82	47.24%
MAISON	84,528.00	68,760.30	52.76%

3.1.6 `NB_PIECES` :

variable numérique indiquant le nombre de pièces du logement

NB_PIECES	nb_lignes	EXPO	pct_EXPO
0	3,718.00	2,413.72	1.85%
1	40,300.00	31,943.51	24.51%
2	60,973.00	51,931.99	39.85%
3	31,695.00	27,730.87	21.28%
4	9,615.00	8,495.12	6.52%
NA	9,350.00	7,807.90	5.99%

3.1.7 `SITUATION_JURIDIQUE` :

variable catégorielle indiquant si le souscripteur est prorpiétaire (PROPRIO) ou locataire (LOCATAIRE) du logement assuré

SITUATION_JURIDIQUE	nb_lignes	EXPO	pct_EXPO
LOCATAIRE	86,421.00	76,135.03	58.42%
PROPRIO	69,230.00	54,188.09	41.58%

3.1.8 `NIVEAU_JURIDIQUE` :

variable catégorielle codant le niveau de couverture de la garantie juridique

NIVEAU_JURIDIQUE	nb_lignes	EXPO	pct_EXPO
JUR1	153,568.00	128,854.03	98.87%
JUR2	2,083.00	1,469.09	1.13%

3.1.9 `VALEUR_DES_BIENS` :

variable numérique reflétant la valeur couverte du contenu du logement

VALEUR_DES_BIENS	nb_lignes	EXPO	pct_EXPO
0	8,307.00	5,259.05	4.04%
3500	51,879.00	40,929.10	31.41%
9000	34,770.00	30,059.53	23.07%
20000	33,745.00	29,664.65	22.76%
35000	11,528.00	10,299.00	7.90%
50000	7,944.00	7,243.56	5.56%
80000	6,534.00	6,018.46	4.62%
100000	944.00	849.77	0.65%

3.1.10 `OBJETS_DE_VALEUR` :

variable catégorielle codant le niveau de couverture des objets

OBJETS_DE_VALEUR	nb_lignes	EXPO	pct_EXPO
NIVEAU_1	145,959.00	121,738.34	93.41%
NIVEAU_2	9,692.00	8,584.78	6.59%

3.1.11 `ZONIER` :

variable catégorielle codant la zone du bien assuré. Le code est constitué d’une lettre représentant une zone et d’un nombre représentant une partie de la zone

ZONIER	nb_lignes	EXPO	pct_EXPO
A1	1,822.00	1,579.15	1.21%
A10	1,125.00	963.04	0.74%
A11	4,775.00	4,140.40	3.18%
A12	3,115.00	2,694.44	2.07%
A13	1,879.00	1,593.61	1.22%
A14	828.00	695.12	0.53%
A2	1,750.00	1,466.96	1.13%
A3	1,129.00	975.75	0.75%
A4	81.00	62.13	0.05%
A5	16.00	14.27	0.01%
A6	698.00	610.63	0.47%
A7	256.00	222.42	0.17%
A8	66.00	61.27	0.05%
A9	4,082.00	3,526.65	2.71%
B1	793.00	701.39	0.54%
B10	3,520.00	3,005.74	2.31%
B11	11.00	8.94	0.01%
B12	342.00	289.71	0.22%
B13	44.00	37.32	0.03%
B14	116.00	98.55	0.08%
B16	211.00	185.81	0.14%
B17	3,530.00	2,992.31	2.30%
B18	8.00	6.25	0.00%
B19	950.00	800.99	0.61%
B2	438.00	385.83	0.30%
B20	6.00	5.58	0.00%
B21	44.00	37.82	0.03%
B22	59.00	52.15	0.04%
B23	214.00	183.13	0.14%
B24	135.00	118.09	0.09%
B25	641.00	551.08	0.42%
B26	1,477.00	1,264.44	0.97%
B27	188.00	167.93	0.13%
B28	268.00	228.86	0.18%
B29	4,269.00	3,645.34	2.80%
B3	1,956.00	1,691.09	1.30%
B30	506.00	443.22	0.34%
B31	79.00	71.24	0.05%
B32	7,015.00	5,906.43	4.53%
B33	186.00	167.03	0.13%
B34	708.00	616.63	0.47%
B35	3.00	2.28	0.00%
B36	454.00	377.39	0.29%
B37	150.00	126.09	0.10%
B38	85.00	70.05	0.05%
B39	2,543.00	2,159.86	1.66%
B4	1,868.00	1,603.70	1.23%
B40	5,776.00	4,859.50	3.73%
B41	402.00	341.69	0.26%
B42	29.00	25.69	0.02%
B43	6,038.00	5,066.02	3.89%
B44	6.00	5.09	0.00%
B45	452.00	379.98	0.29%
B5	34.00	30.47	0.02%
B6	122.00	105.72	0.08%
B7	353.00	303.13	0.23%
B8	239.00	214.63	0.16%
B9	608.00	519.74	0.40%
C1	2,760.00	2,325.40	1.78%
C10	899.00	723.85	0.56%
C11	1,322.00	1,067.75	0.82%
C12	2,453.00	1,917.28	1.47%
C13	1,236.00	959.28	0.74%
C14	1,372.00	1,095.67	0.84%
C15	1,202.00	947.23	0.73%
C16	1,272.00	998.16	0.77%
C17	845.00	656.75	0.50%
C18	739.00	605.29	0.46%
C19	2,897.00	2,317.49	1.78%
C2	6,067.00	5,172.98	3.97%
C20	14,581.00	12,217.25	9.37%
C21	4,816.00	4,002.97	3.07%
C22	2,555.00	2,050.02	1.57%
C23	6,214.00	5,119.00	3.93%
C24	1,972.00	1,621.60	1.24%
C3	3,041.00	2,571.18	1.97%
C4	1,713.00	1,444.03	1.11%
C5	6,833.00	5,665.22	4.35%
C6	5,176.00	4,382.99	3.36%
C7	2,103.00	1,736.06	1.33%
C8	5,550.00	4,455.87	3.42%
C9	9,535.00	7,810.05	5.99%

3.1.12 `NBSIN_TYPE1_AN1` :

variable numérique indiquant le nombre de sinistre de type 1 l’année précédente

NBSIN_TYPE1_AN1	nb_lignes	EXPO	pct_EXPO
0	145,434.00	121,118.77	92.94%
1	9,378.00	8,465.38	6.50%
2	765.00	676.86	0.52%
3	74.00	62.11	0.05%

3.1.13 `NBSIN_TYPE1_AN2` :

variable numérique indiquant le nombre de sinistre de type 1 il y a 2 ans

ATTENTION : comme la table ci-dessous le démontre, cette variable présente un problème puisque toutes les lignes ont au moins un sinistre

Comme discuté dans le hackathon on ignore cette variable

NBSIN_TYPE1_AN2	nb_lignes	EXPO	pct_EXPO
1	146,676.00	122,175.11	93.75%
2	8,225.00	7,462.01	5.73%
3	750.00	685.99	0.53%

3.1.14 `NBSIN_TYPE1_AN3` :

variable numérique indiquant le nombre de sinistre de type 1 il y a 3 ans

NBSIN_TYPE1_AN3	nb_lignes	EXPO	pct_EXPO
0	147,569.00	122,889.56	94.30%
1	7,403.00	6,807.03	5.22%
2	625.00	576.62	0.44%
3	54.00	49.91	0.04%

3.1.15 `NBSIN_TYPE2_AN1` :

variable numérique indiquant le nombre de sinistre de type 2 l’année précédente

NBSIN_TYPE2_AN1	nb_lignes	EXPO	pct_EXPO
0	152,587.00	127,610.08	97.92%
1	2,870.00	2,540.14	1.95%
2	180.00	159.97	0.12%
3	14.00	12.93	0.01%

3.1.16 `NBSIN_TYPE2_AN2` :

variable numérique indiquant le nombre de sinistre de type 2 il y a 2 ans

NBSIN_TYPE2_AN2	nb_lignes	EXPO	pct_EXPO
0	136,330.00	114,011.30	87.48%
1	2,058.00	1,861.63	1.43%
2	113.00	103.72	0.08%
3	8.00	8.00	0.01%
NA	17,142.00	14,338.46	11.00%

3.1.17 `NBSIN_TYPE2_AN3` :

variable numérique indiquant le nombre de sinistre de type 2 il y a 3 ans

NBSIN_TYPE2_AN3	nb_lignes	EXPO	pct_EXPO
0	153,501.00	128,360.07	98.49%
1	2,019.00	1,844.33	1.42%
2	125.00	113.36	0.09%
3	6.00	5.37	0.00%

3.1.18 `id` :

identifiant du risque, clef de jointure

3.1.19 `ANNEE` :

variable catégorielle indiquant l’année d’observation du risque : 2016, 2017, 2018

ANNEE	nb_lignes	EXPO	pct_EXPO
2016	54,283.00	45,171.85	34.66%
2017	52,053.00	43,654.97	33.50%
2018	49,315.00	41,496.30	31.84%

3.2 sin_train

Le fichier expo_train.csv contient la liste des contrats en risques historiquement suivis :

sin_train = read.table(file = here('data', 'raw', 'sin_train.csv'), header=T, sep=';', dec=',', encoding = 'UTF-8', stringsAsFactors = F)

datatable(head(sin_train))

str(sin_train)

## 'data.frame':    2693 obs. of  4 variables:
##  $ id   : int  15 277 643 668 730 1816 1852 2061 2270 2322 ...
##  $ NB   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ COUT : num  521.4 3000.2 26.3 462.4 640.9 ...
##  $ ANNEE: int  2017 2016 2018 2017 2018 2018 2016 2018 2018 2017 ...

id : identifiant du risque, clef de jointure avec le fichier expo_train
NB : Nombre de sinistres de la typologie modélisée (inconnue) sur la période d’observation

NB	nb_lignes
1	2,613.00
2	79.00
3	1.00

3.3 fichier combiné

On peut désormais combine les information d’exposition et de sinistres :

mrh <- expo_train %>%
      left_join(sin_train, by =c('id','ANNEE')) %>%
      replace_na(list('NB'=0, 'COUT'=0))

N’ayant pas accèsau fichier de test utilisé lors du KAckathon, on dissocie ici le fichier mrh en un fichier train (70%), un fichier val (20%) et un fichier test (10%)

set.seed(42)

trainIndex<- createDataPartition(mrh$NB>=0, p=.7, list=FALSE, time=1)

mrh_train<-mrh[trainIndex,]
mrh_test<-mrh[ -trainIndex,]

valIndex<- createDataPartition(mrh_test$NB>=0, p=.66, list=FALSE, time=1)
mrh_val<-mrh_test[ valIndex,]
mrh_test<-mrh_test[ -valIndex,]

3.3.1 Taux de sinistres moyen :

tx <- sum(mrh_train$NB)/sum(mrh_train$EXPO)

print(paste0('taux de sinistre moyen annuel : ', round(tx * 100,2), '%'))

## [1] "taux de sinistre moyen annuel : 2.12%"

3.3.2 Taux de sinistre par variable :

3.3.2.1 FORMULE

FORMULE	Freq	nb_lignes	EXPO	pct_EXPO
ALL_INCLUDE	0.57%	4,695.00	3,355.67	3.67%
CONFORT	2.11%	61,372.00	51,858.85	56.79%
ESSENTIEL	2.22%	15,555.00	13,597.27	14.89%
MEDIUM	1.69%	27,334.00	22,501.75	24.64%

ENSEIGNEMENT : Le comportement observé est stable dans le temps. Les niveaux CONFORT et ESSENTIEL semblent avoir un niveau de risque de fréquence équivalent

3.3.2.2 TYPE_RESIDENCE

TYPE_RESIDENCE	Freq	nb_lignes	EXPO	pct_EXPO
PRINCIPALE	2.15%	91,410.00	78,937.92	86.45%
SECONDAIRE	0.75%	17,546.00	12,375.63	13.55%

ENSEIGNEMENT : Le comportement observé est stable dans le temps. Les résidences SECONDAIRE ont systématiquement une fréquence moindre

3.3.2.3 TYPE_HABITATION

TYPE_HABITATION	Freq	nb_lignes	EXPO	pct_EXPO
APPARTEMENT	1.98%	49,749.00	43,093.59	47.19%
MAISON	1.95%	59,207.00	48,219.96	52.81%

ENSEIGNEMENT : Le comportement observé N’est PAS stable dans le temps. Le type d’habitation est à exclure.

3.3.2.4 NB_PIECES

NB_PIECES	Freq	nb_lignes	EXPO	pct_EXPO
0	0.12%	2,610.00	1,715.52	1.88%
1	1.61%	28,256.00	22,425.90	24.56%
2	1.83%	42,724.00	36,390.78	39.85%
3	2.48%	22,181.00	19,440.34	21.29%
4	2.95%	6,639.00	5,879.31	6.44%
NA	1.98%	6,546.00	5,461.69	5.98%

ENSEIGNEMENT : Le NB_PIECES présente des valeurs manquantes. La fréquence tend à augmenter avec le nombre de pièces.

3.3.2.5 SITUATION_JURIDIQUE

SITUATION_JURIDIQUE	Freq	nb_lignes	EXPO	pct_EXPO
LOCATAIRE	2.21%	60,339.00	53,225.98	58.29%
PROPRIO	1.63%	48,617.00	38,087.57	41.71%

ENSEIGNEMENT : Les propriétaires semblent avoir un fréquence moindre que les locataires de manière consistante, mais avec de la variabilité.

3.3.2.6 NIVEAU_JURIDIQUE

NIVEAU_JURIDIQUE	Freq	nb_lignes	EXPO	pct_EXPO
JUR1	1.97%	107,497.00	90,273.20	98.86%
JUR2	1.52%	1,459.00	1,040.35	1.14%

ENSEIGNEMENT : Le NIVEAU_JURIDIQUE N’est PAS stable dans le temps. A exclure.

3.3.2.7 VALEUR_DES_BIENS

VALEUR_DES_BIENS	Freq	nb_lignes	EXPO	pct_EXPO
0	0.25%	5,782.00	3,682.96	4.03%
3500	1.36%	36,543.00	28,888.32	31.64%
9000	1.95%	24,283.00	21,014.51	23.01%
20000	2.25%	23,639.00	20,768.83	22.74%
35000	2.86%	8,047.00	7,187.60	7.87%
50000	3.14%	5,470.00	4,991.12	5.47%
80000	3.09%	4,557.00	4,210.03	4.61%
100000	3.72%	635.00	570.18	0.62%

ENSEIGNEMENT : La fréquence augmente globalement avec la valeur des biens selon une relation non linéaire

3.3.2.8 OBJETS_DE_VALEUR

OBJETS_DE_VALEUR	Freq	nb_lignes	EXPO	pct_EXPO
NIVEAU_1	1.84%	102,196.00	85,308.12	93.42%
NIVEAU_2	3.77%	6,760.00	6,005.44	6.58%

ENSEIGNEMENT : Le NIVEAU_2 présente une fréquence significativement plus élevée.

3.3.2.9 ZONIER

ZONIER	Freq	nb_lignes	EXPO	pct_EXPO
A1	1.75%	1,271.00	1,098.88	1.20%
A10	1.21%	768.00	659.36	0.72%
A11	0.62%	3,349.00	2,901.64	3.18%
A12	1.10%	2,182.00	1,896.53	2.08%
A13	0.97%	1,326.00	1,112.93	1.22%
A14	1.81%	596.00	497.59	0.54%
A2	1.14%	1,199.00	1,004.07	1.10%
A3	0.57%	814.00	705.81	0.77%
A4	0.00%	55.00	40.92	0.04%
A5	0.00%	11.00	10.24	0.01%
A6	0.72%	495.00	427.35	0.47%
A7	0.69%	170.00	145.53	0.16%
A8	2.29%	48.00	43.73	0.05%
A9	0.56%	2,885.00	2,488.88	2.73%
B1	0.41%	556.00	492.90	0.54%
B10	1.72%	2,492.00	2,140.03	2.34%
B11	0.00%	6.00	4.80	0.01%
B12	2.37%	249.00	210.74	0.23%
B13	6.23%	37.00	32.09	0.04%
B14	0.00%	83.00	69.97	0.08%
B16	0.76%	147.00	130.79	0.14%
B17	1.79%	2,442.00	2,071.38	2.27%
B18	0.00%	6.00	4.25	0.00%
B19	2.42%	672.00	572.48	0.63%
B2	1.09%	308.00	274.89	0.30%
B20	0.00%	6.00	5.58	0.01%
B21	0.00%	30.00	25.92	0.03%
B22	0.00%	41.00	36.37	0.04%
B23	0.78%	150.00	128.79	0.14%
B24	2.56%	91.00	78.21	0.09%
B25	1.33%	436.00	375.81	0.41%
B26	1.24%	1,031.00	884.64	0.97%
B27	0.84%	133.00	119.69	0.13%
B28	1.35%	173.00	148.10	0.16%
B29	1.83%	2,979.00	2,541.22	2.78%
B3	0.75%	1,394.00	1,210.68	1.33%
B30	0.96%	355.00	310.94	0.34%
B31	0.00%	56.00	51.74	0.06%
B32	1.81%	4,913.00	4,134.65	4.53%
B33	0.85%	134.00	117.90	0.13%
B34	1.23%	518.00	446.21	0.49%
B35	0.00%	3.00	2.28	0.00%
B36	1.33%	319.00	266.18	0.29%
B37	0.00%	108.00	94.32	0.10%
B38	5.38%	67.00	55.78	0.06%
B39	1.37%	1,804.00	1,542.57	1.69%
B4	0.73%	1,324.00	1,134.45	1.24%
B40	2.41%	4,045.00	3,427.84	3.75%
B41	2.79%	282.00	242.19	0.27%
B42	0.00%	22.00	20.55	0.02%
B43	1.72%	4,229.00	3,553.62	3.89%
B44	0.00%	5.00	4.09	0.00%
B45	0.83%	317.00	268.14	0.29%
B5	0.00%	24.00	20.47	0.02%
B6	1.30%	91.00	77.15	0.08%
B7	2.84%	243.00	211.00	0.23%
B8	0.00%	152.00	136.97	0.15%
B9	1.93%	437.00	374.18	0.41%
C1	2.68%	1,923.00	1,619.36	1.77%
C10	2.92%	638.00	515.30	0.56%
C11	3.19%	906.00	739.61	0.81%
C12	3.71%	1,700.00	1,335.15	1.46%
C13	2.09%	882.00	681.10	0.75%
C14	2.77%	962.00	769.95	0.84%
C15	3.70%	851.00	666.75	0.73%
C16	2.56%	894.00	707.78	0.78%
C17	2.18%	592.00	463.76	0.51%
C18	1.85%	526.00	429.41	0.47%
C19	3.29%	2,006.00	1,607.93	1.76%
C2	2.34%	4,216.00	3,597.31	3.94%
C20	3.06%	10,231.00	8,576.45	9.39%
C21	1.62%	3,365.00	2,789.43	3.05%
C22	1.50%	1,789.00	1,430.49	1.57%
C23	2.15%	4,418.00	3,642.23	3.99%
C24	1.91%	1,360.00	1,122.90	1.23%
C3	1.18%	2,127.00	1,817.00	1.99%
C4	3.82%	1,174.00	974.79	1.07%
C5	1.79%	4,772.00	3,952.35	4.33%
C6	1.59%	3,614.00	3,073.13	3.37%
C7	2.10%	1,442.00	1,195.34	1.31%
C8	1.65%	3,845.00	3,081.29	3.37%
C9	2.93%	6,644.00	5,438.74	5.96%

ENSEIGNEMENT : Le ZONIER présente une fréquence avec beaucoup de variabilité, mais il semble y avoir une croissance par région qu’on peut tester :

REGION	Freq	nb_lignes	EXPO	pct_EXPO
A	0.92%	15,169.00	13,033.47	14.27%
B	1.66%	32,910.00	28,052.53	30.72%
C	2.40%	60,877.00	50,227.55	55.01%

3.3.2.10 NBSIN_TYPE1_AN1

NBSIN_TYPE1_AN1	Freq	nb_lignes	EXPO	pct_EXPO
0	1.83%	101,793.00	84,856.05	92.93%
1	3.70%	6,581.00	5,945.96	6.51%
2	4.04%	530.00	469.28	0.51%
3	9.46%	52.00	42.26	0.05%

ENSEIGNEMENT : La fréquence augmente globalement avec le nombre de sinistres de type 1 l’année précédente. Toutefois, il semble surtout q’il existe un effet d’une indicatrice a eu sinistre ou non.

3.3.2.11 NBSIN_TYPE1_AN3

NBSIN_TYPE1_AN3	Freq	nb_lignes	EXPO	pct_EXPO
0	1.86%	103,309.00	86,105.62	94.30%
1	3.42%	5,171.00	4,767.33	5.22%
2	6.13%	435.00	403.13	0.44%
3	5.34%	41.00	37.46	0.04%

ENSEIGNEMENT : La fréquence augmente globalement avec le nombre de sinistres de type 1 il y a 3 ans. Toutefois, il semble surtout q’il existe un effet d’une indicatrice a eu sinistre ou non.

3.3.2.12 NBSIN_TYPE2_AN1

NBSIN_TYPE2_AN1	Freq	nb_lignes	EXPO	pct_EXPO
0	1.95%	106,796.00	89,403.29	97.91%
1	2.61%	2,034.00	1,797.95	1.97%
2	2.93%	115.00	102.38	0.11%
3	0.00%	11.00	9.93	0.01%

ENSEIGNEMENT : Il n’existe pas de lien stable sur la fréquence. A ne pas retenir.

3.3.2.13 NBSIN_TYPE2_AN2

NBSIN_TYPE2_AN2	Freq	nb_lignes	EXPO	pct_EXPO
0	1.96%	95,382.00	79,861.28	87.46%
1	3.02%	1,452.00	1,307.09	1.43%
2	1.39%	77.00	71.91	0.08%
3	0.00%	3.00	3.00	0.00%
NA	1.86%	12,042.00	10,070.27	11.03%

ENSEIGNEMENT : Il n’existe pas de lien stable sur la fréquence. A ne pas retenir.

3.3.2.14 NBSIN_TYPE2_AN3

NBSIN_TYPE2_AN3	Freq	nb_lignes	EXPO	pct_EXPO
0	1.95%	107,445.00	89,931.57	98.49%
1	2.68%	1,414.00	1,292.12	1.42%
2	4.31%	93.00	85.86	0.09%
3	0.00%	4.00	4.00	0.00%

ENSEIGNEMENT : Il n’existe pas de lien stable sur la fréquence. A ne pas retenir.

3.3.2.15 ANNEE

ANNEE	Freq	nb_lignes	EXPO	pct_EXPO
2016	1.82%	37,947.00	31,635.56	34.64%
2017	2.04%	36,440.00	30,597.37	33.51%
2018	2.04%	34,569.00	29,080.62	31.85%

ENSEIGNEMENT : L’année 2016 semble atypique. Le niveau actuel semble plus proche de ceux de 2017 et 2018.

3.3.3 Nettoyage des données

preprocess<- function(dt) {
            res <- dt %>%
            mutate( 
              NB_PIECES = ifelse(is.na(NB_PIECES), 2, NB_PIECES), #2 est le mode de NB_PIECES
              NBSIN_TYPE1_AN1 = ifelse(NBSIN_TYPE1_AN1==0,0,1), #On regroupe les modalités non stables dans le temps
              NBSIN_TYPE1_AN3 = ifelse(NBSIN_TYPE1_AN3==0,0,1), #On regroupe les modalités non stables dans le temps
              REGION =  gsub('[0-9]', '', ZONIER) #création de région
            ) %>%
            select(-c('X', 'TYPE_HABITATION', 'NIVEAU_JURIDIQUE', 'NBSIN_TYPE1_AN2', 'NBSIN_TYPE2_AN1', 'NBSIN_TYPE2_AN2', 'NBSIN_TYPE2_AN3'  )) %>% # On supprime les variables sans lien stable dans le temps
            relocate('id', 'EXPO', 'FORMULE', 'TYPE_RESIDENCE', 'SITUATION_JURIDIQUE', 'OBJETS_DE_VALEUR', 'VALEUR_DES_BIENS', 'NB_PIECES', 'ZONIER', 'REGION', 'NBSIN_TYPE1_AN1', 'NBSIN_TYPE1_AN3', 'ANNEE', 'NB', 'COUT' )

            return(res)
}

df_train <- preprocess(mrh_train)
df_val <- preprocess(mrh_val)
df_test <- preprocess(mrh_test)

3.3.4 Calcul du V de Cramer

On mesude l’association des variables entre elles

X_train <- df_train %>% select(-c('id', 'EXPO', 'NB', 'COUT'))


# Initialize empty matrix to store coefficients
empty_m <- matrix(ncol = length(X_train),
            nrow = length(X_train),
            dimnames = list(names(X_train), 
                            names(X_train)))
# Function that accepts matrix for coefficients and data and returns a correlation matrix
calculate_cramer <- function(m, df) {
 for (r in seq(nrow(m))){
   for (c in seq(ncol(m))){
     m[[r, c]] <- lsr::cramersV(X_train[[r]], X_train[[c]])
   }
 }
    return(m)
}

cor_matrix <- calculate_cramer(empty_m ,X_train)

corrplot(cor_matrix)

4 Modélisation de la fréquence

Dans cette partie, on essaye de modéliser la fréquence des sinistres.

On commence par une modélisation de référence qui consistera en un modèle GLM Poisson classique en tirant partie des enseignements de la partie précédente.

h2o.init(nthreads = -1)

##  Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         1 hours 46 minutes 
##     H2O cluster timezone:       Europe/Paris 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.32.1.3 
##     H2O cluster version age:    1 month and 28 days  
##     H2O cluster name:           H2O_started_from_R_fabien_yji981 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   29.89 GB 
##     H2O cluster total cores:    24 
##     H2O cluster allowed cores:  24 
##     H2O cluster healthy:        TRUE 
##     H2O Connection ip:          localhost 
##     H2O Connection port:        54321 
##     H2O Connection proxy:       NA 
##     H2O Internal Security:      FALSE 
##     H2O API Extensions:         Amazon S3, XGBoost, Algos, AutoML, Core V3, TargetEncoder, Core V4 
##     R Version:                  R version 4.1.0 (2021-05-18)

reg_bst_100 = h2o.gbm(y = 'NB', x = names(df_train),
                       distribution = "poisson",
                       offset_column = "EXPO",
                       training_frame = as.h2o(df_train),
                       validation_frame = as.h2o(df_val),
                       ntrees = 100,
                       nfolds = 5,
                       seed = 1)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================| 100%

plot(reg_bst_100)

reg0 = glm(NB~1+offset(log(EXPO)),family=poisson,data=df_train)

summary(reg0)

## 
## Call:
## glm(formula = NB ~ 1 + offset(log(EXPO)), family = poisson, data = df_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.2058  -0.2058  -0.2058  -0.1774   4.8750  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.85523    0.02274  -169.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 15507  on 108955  degrees of freedom
## Residual deviance: 15507  on 108955  degrees of freedom
## AIC: 19294
## 
## Number of Fisher Scoring iterations: 6

plot(reg0)

fit_glm <- glm(NB ~ FORMULE + TYPE_RESIDENCE + SITUATION_JURIDIQUE + OBJETS_DE_VALEUR + as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) + REGION + as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) + as.factor(ANNEE), 
              data=df_train, offset=log(EXPO), family=quasipoisson())

summary(fit_glm)

## 
## Call:
## glm(formula = NB ~ FORMULE + TYPE_RESIDENCE + SITUATION_JURIDIQUE + 
##     OBJETS_DE_VALEUR + as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) + 
##     REGION + as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) + 
##     as.factor(ANNEE), family = quasipoisson(), data = df_train, 
##     offset = log(EXPO))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4944  -0.2199  -0.1767  -0.1221   4.3198  
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -8.35092    0.78731 -10.607  < 2e-16 ***
## FORMULECONFORT                    0.89236    0.21816   4.090 4.31e-05 ***
## FORMULEESSENTIEL                  0.85689    0.22565   3.798 0.000146 ***
## FORMULEMEDIUM                     0.85491    0.21928   3.899 9.67e-05 ***
## TYPE_RESIDENCESECONDAIRE         -0.61505    0.11032  -5.575 2.48e-08 ***
## SITUATION_JURIDIQUEPROPRIO       -0.05789    0.06008  -0.963 0.335315    
## OBJETS_DE_VALEURNIVEAU_2          0.28047    0.08372   3.350 0.000808 ***
## as.factor(VALEUR_DES_BIENS)3500   0.62766    0.34944   1.796 0.072467 .  
## as.factor(VALEUR_DES_BIENS)9000   0.77136    0.35190   2.192 0.028381 *  
## as.factor(VALEUR_DES_BIENS)20000  0.83321    0.35253   2.364 0.018103 *  
## as.factor(VALEUR_DES_BIENS)35000  0.97408    0.35808   2.720 0.006523 ** 
## as.factor(VALEUR_DES_BIENS)50000  0.98658    0.36156   2.729 0.006360 ** 
## as.factor(VALEUR_DES_BIENS)80000  0.84862    0.36617   2.318 0.020477 *  
## as.factor(VALEUR_DES_BIENS)1e+05  0.92804    0.42120   2.203 0.027574 *  
## as.factor(NB_PIECES)1             2.04091    0.82752   2.466 0.013653 *  
## as.factor(NB_PIECES)2             1.98259    0.82633   2.399 0.016430 *  
## as.factor(NB_PIECES)3             2.23114    0.82809   2.694 0.007054 ** 
## as.factor(NB_PIECES)4             2.32381    0.83140   2.795 0.005190 ** 
## REGIONB                           0.43795    0.10544   4.154 3.27e-05 ***
## REGIONC                           0.95079    0.09975   9.532  < 2e-16 ***
## as.factor(NBSIN_TYPE1_AN1)1       0.43828    0.07264   6.034 1.61e-09 ***
## as.factor(NBSIN_TYPE1_AN3)1       0.33269    0.08157   4.078 4.54e-05 ***
## as.factor(ANNEE)2017              0.10300    0.05889   1.749 0.080290 .  
## as.factor(ANNEE)2018              0.10260    0.05960   1.721 0.085171 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.119526)
## 
##     Null deviance: 15507  on 108955  degrees of freedom
## Residual deviance: 14905  on 108932  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 9

plot(fit_glm)

x <- c('FORMULE', 'TYPE_RESIDENCE', 'SITUATION_JURIDIQUE','OBJETS_DE_VALEUR', 'VALEUR_DES_BIENS', 'NB_PIECES', 'REGION', 'NBSIN_TYPE1_AN1', 'NBSIN_TYPE1_AN3', 'ANNEE')
y <- 'NB'
w <- 'EXPO'

# Input maker
prep_xgb <- function(dat, x) {
  data.matrix(dat[, x, drop = FALSE])
}
# Data interface to XGBoost
dtrain <- xgb.DMatrix(
  prep_xgb(df_train, x), 
  label = df_train[[y]], 
  weight = df_train[[w]]
)
# Parameters chosen by 5-fold grouped CV
params_freq <- list(
  learning_rate = 0.2,
  max_depth = 5,
  alpha = 3,
  lambda = 0.5,
  max_delta_step = 2,
  min_split_loss = 0,
  colsample_bytree = 1,
  subsample = 0.9
)
# Fit
set.seed(1)
fit_xgb <- xgb.train(
  params_freq, 
  data = dtrain,
  nrounds = 580,
  objective = "count:poisson",
  watchlist = list(train = dtrain),
  print_every_n = 100
)

## [1]  train-poisson-nloglik:0.501726 
## [101]    train-poisson-nloglik:0.108273 
## [201]    train-poisson-nloglik:0.093639 
## [301]    train-poisson-nloglik:0.092993 
## [401]    train-poisson-nloglik:0.092593 
## [501]    train-poisson-nloglik:0.092233 
## [580]    train-poisson-nloglik:0.091997

# Save and load model
xgb.save(fit_xgb, "xgb.model")

## [1] TRUE

fit_xgb <- xgb.load("xgb.model")

5 Model Explanations

The models are ready, so let’s shed light into them.

5.1 Setting up explainers

We start by setting up the explainers. These are basically objects that know how to create predictions.

Crucial: the prediction function needs to work for subsets of datasets, not just for the original model data set.

fl_glm <- flashlight(
  model = fit_glm, label = "GLM", 
  predict_function = function(fit, X) predict(fit, X, type = "response")
)

fl_xgb <- flashlight(
  model = fit_xgb, label = "XGBoost", 
  predict_function = function(fit, X) predict(fit, prep_xgb(X, x))
)

fl_xgb <- flashlight(
  model = fit_xgb, label = "H20 GBM", 
  predict_function = function(fit, X) predict(fit, prep_xgb(X, x))
)

reg_bst_100

## Model Details:
## ==============
## 
## H2ORegressionModel: gbm
## Model ID:  GBM_model_R_1626594240864_7 
## Model Summary: 
##   number_of_trees number_of_internal_trees model_size_in_bytes min_depth
## 1             100                      100               27915         5
##   max_depth mean_depth min_leaves max_leaves mean_leaves
## 1         5    5.00000          7         28    17.52000
## 
## 
## H2ORegressionMetrics: gbm
## ** Reported on training data. **
## 
## MSE:  0.0006816525
## RMSE:  0.02610848
## MAE:  0.002147957
## RMSLE:  0.01218481
## Mean Residual Deviance :  0.03465347
## 
## 
## H2ORegressionMetrics: gbm
## ** Reported on validation data. **
## 
## MSE:  0.001115619
## RMSE:  0.03340088
## MAE:  0.002735446
## RMSLE:  0.01622289
## Mean Residual Deviance :  0.03836132
## 
## 
## H2ORegressionMetrics: gbm
## ** Reported on cross-validation data. **
## ** 5-fold cross-validation on training data (Metrics computed for combined holdout predictions) **
## 
## MSE:  0.000999386
## RMSE:  0.03161307
## MAE:  0.00255466
## RMSLE:  0.01495924
## Mean Residual Deviance :  0.03593166
## 
## 
## Cross-Validation Metrics Summary: 
##                                mean           sd   cv_1_valid   cv_2_valid
## mae                    0.0025544818 8.2643186E-5  0.002471777 0.0026899083
## mean_residual_deviance   0.03593359 0.0025269587   0.03283671  0.034950383
## mse                     9.992188E-4  8.343817E-5 0.0010756224 0.0011041489
## r2                       0.94598985 0.0046124477    0.9399748    0.9421698
## residual_deviance        0.03593359 0.0025269587   0.03283671  0.034950383
## rmse                    0.031588703 0.0013098742   0.03279668  0.033228736
## rmsle                   0.014955509 3.3117784E-4  0.014847655  0.015446129
##                          cv_3_valid   cv_4_valid   cv_5_valid
## mae                    0.0025595813 0.0025404296 0.0025107136
## mean_residual_deviance   0.03734176    0.0350979   0.03944121
## mse                    9.3566976E-4 9.3680545E-4  9.438475E-4
## r2                       0.94894713   0.94844747    0.9504102
## residual_deviance        0.03734176    0.0350979   0.03944121
## rmse                     0.03058872  0.030607278  0.030722102
## rmsle                   0.014883362  0.014542394  0.015058007

# Combine them and add common elements like reference data
metrics <- list(`Average deviance` = deviance_poisson, 
                `Relative deviance reduction` = r_squared_poisson)
fls <- multiflashlight(list(fl_glm, fl_xgb), data = df_val, 
                       y = y, w = w, metrics = metrics)
# Version on canonical scale
fls_log <- multiflashlight(fls, linkinv = log)

5.2 Performance

We start to interpret the models by looking at model performance using deviance related metrics.

fillc <- "#E69F00"
perf <- light_performance(fls)
perf

## 
## I am an object of class light_performance_multi 
## 
## data.frames:
## 
##  data 
## # A tibble: 4 x 3
##   label   metric                       value
##   <fct>   <fct>                        <dbl>
## 1 GLM     Average deviance            0.151 
## 2 GLM     Relative deviance reduction 0.0506
## 3 H20 GBM Average deviance            0.152 
## 4 H20 GBM Relative deviance reduction 0.0454

plot(perf, geom = "point") +
  labs(x = element_blank(), y = element_blank())

5.3 Importance

Next, we consider permutation variable importance. By default flashlight uses the first performance metric specified above.

imp <- light_importance(fls, v = x)
plot(imp, fill = fillc, color = "black")

5.3.1 Partial dependence curves

Taking the average of many ICE curves produces the famous partial dependence plot. We consider such plots for four predictors.

plot(light_profile(fls, v = "VALEUR_DES_BIENS"))

plot(light_profile(fls, v = "NB_PIECES"))

plot(light_profile(fls, v = "TYPE_RESIDENCE"))

5.4 Classic diagnostic plots

Classic predicted/residual/response versus covariable plots are worth a look.

# Average predicted versus covariable
plot(light_profile(fls, v = "NB_PIECES", type = "predicted"))

# Average residual versus covariable
plot(light_profile(fls, v = "NB_PIECES", type = "residual")) +
  geom_hline(yintercept = 0)

# Average response versus covariable
plot(light_profile(fls, v = "NB_PIECES", type = "response"))

5.4.1 Multiple aspects combined

We often get a good picture of the effect of a covariable by combining partial dependence with classic plots.

eff_DrivAge <- light_effects(fls, v = "NB_PIECES", counts_weighted = TRUE)
p <- plot(eff_DrivAge, show_points = FALSE)
plot_counts(p, eff_DrivAge, alpha = 0.3)

interact_rel <- light_interaction(
  fls_log, 
  v = most_important(imp, 4), 
  take_sqrt = FALSE,
  pairwise = TRUE, 
  use_linkinv = TRUE,
  seed = 61
)
plot(interact_rel, color = "black", fill = fillc, rotate_x = TRUE)

5.4.2 On absolute scale

interact_abs <- light_interaction(
  fls_log, 
  v = most_important(imp, 4), 
  normalize = FALSE,
  pairwise = TRUE, 
  use_linkinv = TRUE,
  seed = 61
)
plot(interact_abs, color = "black", fill = fillc, rotate_x = TRUE)

5.4.3 Visualization

Let’s illustrate a strong and a weak interaction by conditional partial dependence plots.

# Strong interaction
# p <- plot(eff_DrivAge, show_points = FALSE)
# plot_counts(p, eff_DrivAge, alpha = 0.3)


pdp_NBPIECES_REGION <- light_profile(fls_log, v = "NB_PIECES", by = "REGION", 
                                  pd_seed = 50, data = df_val, counts_weighted=TRUE)
p <- plot(pdp_NBPIECES_REGION, show_points = FALSE)
pdp_count <- light_effects(fls, v = "NB_PIECES")
plot_counts(p, pdp_count, alpha = 0.3)

pdp_VALEURDESBIENS_REGION <- light_profile(fls_log, v = "VALEUR_DES_BIENS", by = "REGION", 
                                  pd_seed = 50, data = df_val, counts_weighted=TRUE)
plot(pdp_VALEURDESBIENS_REGION)

# Weak interaction
pdp_TYPE_RESIDENCE_REGION <- light_profile(fls_log, v = "TYPE_RESIDENCE", 
                                 by = "REGION", pd_seed = 50, data = df_val, counts_weighted=TRUE)
plot(pdp_TYPE_RESIDENCE_REGION)

summary(df_val$NB/df_val$EXPO)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.00000  0.00000  0.00000  0.02384  0.00000 52.28575

summary(fl_glm$predict_function(fit_glm, df_val))

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.590e-06 8.068e-03 1.654e-02 1.781e-02 2.451e-02 1.296e-01

summary(fl_xgb$predict_function(fit_xgb, df_val))

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.001192 0.009000 0.017827 0.019045 0.025425 0.171986

fls

DSA 2021 - TD R

Fabien FAIVRE

18 juillet 2021

1 Présentation du projet

2 Analyse

2.1 Mise en place de l’environnement

2.2 code

3 Présentation du jeu de données

3.1 expo_train

3.1.1 X:

3.1.2 EXPO :

3.1.3 FORMULE :

3.1.4 TYPE_RESIDENCE :

3.1.5 TYPE_HABITATION :

3.1.6 NB_PIECES :

3.1.7 SITUATION_JURIDIQUE :

3.1.8 NIVEAU_JURIDIQUE :

3.1.9 VALEUR_DES_BIENS :

3.1.10 OBJETS_DE_VALEUR :

3.1.11 ZONIER :

3.1.12 NBSIN_TYPE1_AN1 :

3.1.13 NBSIN_TYPE1_AN2 :

3.1.14 NBSIN_TYPE1_AN3 :

3.1.15 NBSIN_TYPE2_AN1 :

3.1.16 NBSIN_TYPE2_AN2 :

3.1.17 NBSIN_TYPE2_AN3 :

3.1.18 id :

3.1.19 ANNEE :

3.2 sin_train

3.3 fichier combiné

3.3.1 Taux de sinistres moyen :

3.3.2 Taux de sinistre par variable :

3.3.2.1 FORMULE

3.3.2.2 TYPE_RESIDENCE

3.3.2.3 TYPE_HABITATION

3.3.2.4 NB_PIECES

3.3.2.5 SITUATION_JURIDIQUE

3.3.2.6 NIVEAU_JURIDIQUE

3.3.2.7 VALEUR_DES_BIENS

3.3.2.8 OBJETS_DE_VALEUR

3.3.2.9 ZONIER

3.3.2.10 NBSIN_TYPE1_AN1

3.3.2.11 NBSIN_TYPE1_AN3

3.3.2.12 NBSIN_TYPE2_AN1

3.3.2.13 NBSIN_TYPE2_AN2

3.3.2.14 NBSIN_TYPE2_AN3

3.3.2.15 ANNEE

3.3.3 Nettoyage des données

3.3.4 Calcul du V de Cramer

4 Modélisation de la fréquence

5 Model Explanations

5.1 Setting up explainers

5.2 Performance

5.3 Importance

5.3.1 Partial dependence curves

5.4 Classic diagnostic plots

5.4.1 Multiple aspects combined

5.4.2 On absolute scale

5.4.3 Visualization

3.1.1 `X`:

3.1.2 `EXPO` :

3.1.3 `FORMULE` :

3.1.4 `TYPE_RESIDENCE` :

3.1.5 `TYPE_HABITATION` :

3.1.6 `NB_PIECES` :

3.1.7 `SITUATION_JURIDIQUE` :

3.1.8 `NIVEAU_JURIDIQUE` :

3.1.9 `VALEUR_DES_BIENS` :

3.1.10 `OBJETS_DE_VALEUR` :

3.1.11 `ZONIER` :

3.1.12 `NBSIN_TYPE1_AN1` :

3.1.13 `NBSIN_TYPE1_AN2` :

3.1.14 `NBSIN_TYPE1_AN3` :

3.1.15 `NBSIN_TYPE2_AN1` :

3.1.16 `NBSIN_TYPE2_AN2` :

3.1.17 `NBSIN_TYPE2_AN3` :

3.1.18 `id` :

3.1.19 `ANNEE` :